from sympy.vector import CoordSys3D, divergence
from sympy import *
import math

R=CoordSys3D('R')
pi=Symbol('Pi')
lam=1
mu=2
u1=sin(pi*R.x)*sin(pi*R.y)
u2=R.x*(R.x-1)*R.y*(R.y-1)
sigma_u1=(lam*diff(u1,R.x,1,R.y,0)+2*mu*diff(u1,R.x,1,R.y,0)+lam*diff(u2,R.x,0,R.y,1))*R.i+(mu*diff(u1,R.x,0,R.y,1)+mu*diff(u2,R.x,1,R.y,0))*R.j
#print(sigma_u1)
#print((divergence(sigma_u1)))
div_sig_u1=divergence(sigma_u1)
f1=-(lam+2*mu)*(-pi**2*sin(pi*R.x)*sin(pi*R.y))-(lam+mu)*((2*R.x-1)*(2*R.y-1))-mu*(-pi**2*sin(pi*R.x)*sin(pi*R.y))

print(div_sig_u1-f1)



